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Abstract  (cont.) 

I  propose  to  ],ook  at  the  simplest  non— trivial  form  of  conservative 
systems,  the  area-preserving  maps,  to  provide  a  new  source  of  examples 
for  investigation  into  the  fundamental  issues  of  descriptive  language, 
style  of  reasoning,  and  representation  techniques.)  1  study  how  physicists 
explore  the  dynamics  of  these  systems'^  jud4e±otfs  use  of  numerical 
experiments.**  To  automate  the  experimenting  process,  two  key  problems  - 
automatic  experiment  control,  and  -f?? “result  interpretation  -  have 
to  be  solved.  Knowledge  of  qualitative  dynamics  and  bifurcations  provides 
a  strong  constraint  on  the  type  of  dynamical  behavior  possible.  This 
constraint  can  be  exploited  to  solve  the  problems .^3-— discusa-  an  approach 
to  the  control  problem Abased  on  this  idea.  The  main  result  paper L 

is  an  implemented  program  which  solves  the  interpretation  problem  by 
using  techniques  from  computational  geometry  and  computer  vision. 


|  Accession  For 

_NTIS  GRA&I 
DTIC  TAB 
Unannounced 
Justif lent!  on _ 


By - 

Distribution/ 

Availability  Codas 
A /ail  and/or 
:Cist  (  Spocial 


MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 
ARTIFICIAL  INTELLIGENCE  LABORATORY 


i 

A. I.  Memo  No.  950 


March  1987 


Extracting  Qualitative  Dynamics  from 
Numerical  Experiments 

Kenneth  Man-kam  Yip 


Abstract 


One  central  problem  in  Qualitative  Physics  is  the  qualitative  prediction  of  long¬ 
time  behavior  of  physical  dynamic  systems.  The  machinery  developed  for  qualitative 
reasoning  -  qualitative  state  vector,  quantity  space,  and  limit  analysis  -  are  largely 
applicable  to  systems  which  are  piecewise  well-approximated  by  low-order  linear  sys¬ 
tems  or  by  first  order  nonlinear  differential  equations.  Typical  nonlinear  systems  - 
those  commonly  encountered  in  Physics  -  exhibit  a  far  richer  spectrum  of  dynamical 
behavior.  I  propose  to  look  at  the  simplest  non-trivial  form  of  conservative  systems, 
the  area- preserving  maps,  to  provide  a  new  source  of  examples  for  investigation  into 
the  fundamental  issues  of  descriptive  language,  style  of  reasoning,  and  representation 
techniques.  I  study  how  physicists  explore  the  dynamics  of  these  systems  by  judi¬ 
cious  use  of  numerical  experiments.  To  automate  the  experimenting  process,  two  key 
problems  -  (1)  automatic  experiment  control,  and  (2)  result  interpretation  -  have  to 
be  solved.  Knowledge  of  qualitative  dynamics  and  bifurcations  provides  a  strong  con¬ 
straint  on  the  type  of  dynamical  behavior  possible.  This  constraint  can  be  exploited 
to  solve  the  problems.  I  discuss  an  approach  to  the  control  problem  based  on  this 
idea.  Finally,  I  describe  the  main  result  of  this  paper:  an  implemented  program  which 
solves  the  interpretation  problem  by  using  techniques  from  computational  geometry 
and  computer  vision. 
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Everything  should  be  made  as  simple  as  possible, 

BUT  NOT  SIMPLER. 
-  A.  Einstein 


1  INTRODUCTION 

Qualitative  Physics  is  a  young  field.  Progress  is  made  when  researchers  formalize 
and  implement  their  understanding  of  how  certain  qualitative  reasoning  tasks,  such 
as  prediction  of  future  behavior,  and  explanation  of  how  the  behavior  comes  about, 
are  being  performed  in  particular  problem  domains.  Two  domains,  among  others, 
have  received  much  attention:  circuit  analysis  and  design  in  the  engineering  domain, 
and  simple  boilers  and  fluid  flow  in  commonsense  physics.  Early  works  in  Qualitative 
Physics  primarily  dealt  with  incremental  deviation  from  equilibrium  states  where  time 
evolution  is  not  explicitly  considered  [4].  More  recent  works  attempt  to  extend  DeK- 
leer’s  qualitative  algebra  and  incremental  analysis  to  handle  time-varying  behavior 
[5,15,14,11]. 

A  common  characteristic  of  the  examples  examined  so  far  is  this:  the  physical 
dynamic  systems  studied  in  these  disciplines  can  often  be  divided  into  different  op¬ 
erating  regions  each  of  which  can  be  modeled  to  a  first  approximation  by  simple  low 
order  linear  differential  equations,  or  by  first  order  nonlinear  differential  equations. 
The  behavior  of  linear  systems  is  particularly  simple:  the  complete  input-output  be¬ 
havior  can  be  summarized  in  a  single  system  transfer  function.  Consequently,  if  the 
response  to  one  type  of  input  is  known,  no  more  information  is  needed  to  determine 
responses  for  other  input  signals. 

The  situation  in  a  nonlinear  system  is  completely  different:  essential  changes  in 
the  qualitative  behavior  of  the  system  may  occur  as  the  amplitude  of  the  input  signal 
changes,  or  as  the  starting  conditions  are  varied.  More  importantly,  nonlinear  systems 
have  a  far  richer  spectrum  of  dynamical  behavior.  Simple  equilibrium  points,  periodic 
and  quasiperiodic  motion,  limit  cycles,  chaotic  motion  as  unpredictable  as  a  sequence 
of  coin  tosses  -  these  are  some  of  the  behavior  found  in  a  typical  nonlinear  system. 

Unfortunately,  these  nonlinear  characteristics  do  not  show  up  in  first  order  non¬ 
linear  differential  equations.  This  is  because  the  continuity  and  (local)  uniqueness  of 
flow  severely  constrain  the  kind  of  behavior  possible  on  the  real  line:  the  flow  either 
tends  towards  an  equilibrium,  or  goes  off  to  infinity. 


In  this  research,  I  therefore  propose  to  look  at  dynamical  systems  -  those  typically 
encountered  in  Physics  -  to  provide  a  new  source  of  examples  for  investigation  into 
the  fundamental  issues  of  descriptive  language,  style  of  reasoning,  and  representation 
techniques  in  qualitative  reasoning  about  nonlinear  dynamical  systems.  Specifically,  1 
will  consider  two-dimensional  discrete  dynamical  systems  defined  by  area- preserving 
maps  containing  a  single  control  parameter.  The  study  of  area-preserving  maps  - 
transformations  of  the  plane  which  preserves  area  -  began  with  the  venerable  prob¬ 
lem  of  the  stability  of  the  solar  system.  I  choose  to  investigate  this  simplest  non-trivial 
type  of  conservative  system  because  many  important  problems  in  physics  -  the  re¬ 
stricted  3-body  problem,  orbits  of  particles  in  accelerators,  and  two  coupled  nonlinear 
oscillators,  just  to  mention  a  few  -  can  be  reduced  to  the  study  of  area-preserving 
maps. 

To  illustrate  the  vocabularies  that  a  physicist  uses  to  describe  the  qualitative 
behavior  of  nonlinear  systems,  we  can  examine  the  following  description  of  one  of  the 
phase  portraits  obtained  from  the  results  of  numerical  experiments  with  the  quadratic 
map  (see  next  section): 


Figure  la  represents  a  number  of  trajectories  (sequences  of  points)  for  coea  = 
0.4  . . .  a  regular  structure  in  a  neighborhood  of  the  elliptic  fixed  point  at  the 
origin,  and  farther  away  a  chaotic  zone.  In  many  places  the  plotted  points  are 
so  dense  that  they  give  the  illusion  of  a  continuous  curve.  Near  the  origin, 
the  “ curves ”  are  almost  circular.  As  we  move  outward,  the  curves  become 
distorted.  Just  inside  the  outermost  regular  curve  lies  a  chain  of  six  closed 
curves,  or  “chain  of  islands”.  Successive  points  of  a  trajectory  jump  from  one 
island  to  another  by  application  of  the  mapping.  Finally,  as  the  curves  become 
more  distorted,  there  is  a  sudden  break-up  and  the  set  of  points  no  longer  lies 
on  a  curve,  but  seems  rather  to  fill  a  two-dimensional  region. 


The  description  is  an  edited  form  of  a  passage  from  [8,  pages  99-100J.  What  is 
striking  in  this  description  is  that  the  language  is  entirely  geometric:  it  contains  terms 
like  fixed  point,  continuous  curve,  chain  of  islands,  two-dimensional  region  etc.  The 
phase  portrait  is  completely  characterized  by  the  spatial  relationships  among  these 
geometric  objects. 

A  goal  of  this  research  is  to  develop  a  clear  understanding  of  the  geometric  lan¬ 
guage  displayed  in  the  above  description.  In  particular,  I  want  to  understand  how 
representation  and  reasoning  techniques  based  on  the  geometric  language  can  be  used 
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as  a  basis  for  a  program  to  perform  automatic  numerical  experiment  control  and  in¬ 
terpretation. 

The  paper  is  organized  as  follows.  Section  two  defines  the  task.  Section  three 
reviews  the  knowledge  -  qualitative  dynamics  and  bifurcations  -  required  for  the 
task.  Section  four  describes  my  approach  to  the  experiment  control  problem.  This 
part  has  not  been  implemented.  Section  five,  the  main  result  of  this  paper,  shows  the 
major  pieces  of  an  implemented  program  which  solves  the  interpretation  problem. 


2  THE  TASK 


Given  an  one-parameter  area- preserving  map  defining  a  discrete  dynamical  system, 
I  am  interested  in  describing  the  qualitative  long-term,  or  asymptotic,  behavior  of 
the  system  over  its  entire  operating  range.  Specifically,  if  U  is  an  open  subset  of  the 
phase  space  X  C  S  x  S,  and  J  C  3?  an  interval  over  which  the  control  parameter 
varies,  I  want  my  program  to  automatically  generate  (1)  a  family  of  phase  portraits 
describing  the  main  dynamical  properties  of  the  map  for  all  initial  conditions  in  U 
and  parameter  values  in  J,  and  (2)  a  summary  graph  which  shows  a  cross-section  of 
the  map  along  certain  axis  of  symmetry  as  a  function  of  the  parameter. 

To  explore  the  complete  dynamics  of  a  nonlinear  system  over  a  large  region  of  the 
phase  space  and  parameter  space  is  a  fairly  typical  problem  in  the  physics  literature. 
A  good  illustration  of  this  task  is  provided  by  Henon’s  well-known  paper,  “Numerical 
Study  of  Quadratic  Area-Preserving  Mappings”  [9].  The  goal  of  Henon's  paper  is  to 
provide  a  description  of  the  main  properties  of  the  quadratic  map : 

xn+1  =  xntosa  -  (yn  -  x*)sina 

yn+i  =  xnsina  -F  (yn  -  x£)cosa  (1) 

where  x  and  y  are  the  state  variables,  and  a  is  the  control  parameter.  The  main 
results  of  Henon’s  paper  are  shown  in  Figures  l(a)-(f),  which  display  the  output  of 
many  numerical  simulations,  and  Figure  2  which  summarizes  the  dynamics  of  Ihe 
map  in  a  two-dimensional  plot.  Dots  on  the  plot  indicate  regions  of  invariant  curve, 
circles  indicate  islands,  and  blanks  indicate  chaotic  zones. 

The  simplest  approach  to  this  problem  is  the  brute  force  method:  it  divides  the 
phase  space  and  parameter  space  into  small  grids  and  tries  every  possible  combi¬ 
nations  of  initial  conditions  and  parameter  values.  A  simple  calculation  will  show 
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Figure  1:  A  partial  list  of  phase  portraits  from  numerical  experiments,  (a)  a  =  1.16 
(b)  a  =  1.33  (c)  a  =  1.58  (d)  a  =  2.0  (e)  a  =  2.04  (f)  a  =  2.21.  Dashed  line:  axis  of 
symmetry. 

that  this  involves  an  enormous  amount  of  computation.  For  instance,  if  we  choose 
a  uniform  grid  size  of  0.01,  we  have  to  compute  approximately  300  x  300  x  600  = 
54  million  orbits.  Assuming,  on  the  average,  0.02  second  is  needed  to  compute  a 
trajectory  of  500  points,  it  will  take  over  300  hours  of  computation  time  to  compute 
all  the  trajectories. 

The  brute  force  method  suffers  from  two  serious  problems.  First,  it  is  grossly 
inefficient  because  most  of  the  phase  portraits  computed  will  be  qualitatively  the 
same.  Second,  it  is  not  reliable  because  there  is  always  the  danger  of  missing  some 
important  qualitative  features  when  the  change  occurs  at  a  resolution  finer  than  the 
grid  size. 
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Figure  2:  Summary  picture  for  all  values  of  a,  and  initial  points  on  the  axis  of 
symmetry.  Dots:  regions  of  invariant  curve.  Circles:  islands.  Blanks:  chaotic  zones. 
Dashed  curves:  unstable  invariant  points. 

A  physicist  often  does  much  better  than  this.  Figure  3  represents  a  flow-chart  of 
what  an  expert  physicist  does  during  the  numerical  experiment.  The  flow-chart  has 
two  nested  loops.  The  outer  loop  involves  deciding  when  to  stop  the  experiment;  the 
inner  loop,  when  to  move  on  to  next  parameter  value.  Controlling  what  experiment 
to  do  next,  and  interpreting  the  results  of  the  simulation  -  these  are  the  two  most 
important  decisions  the  experimenter  has  to  make. 

The  task  of  behavior  prediction  can  now  be  summarized  as  this:  to  develop  a 
picture  of  all  possible  solutions  to  the  dynamical  system  from  a  limited  amount  of 
numerical  experiments  at  a  limited  number  of  initial  conditions  and  parameter  values. 
The  key  observation  is  that  knowledge  of  qualitative  dynamics  and  their  geometric 
manifestations  in  the  phase  space  provides  a  strong  constraint  on  the  type  of  behavior 
possible.  As  we  will  see  in  the  next  section,  this  constraint  translates  into  a  dramatic 
reduction  of  the  amount  of  search  required  to  find  those  combinations  of  states  and 
parameter  values  that  lead  to  “interesting”  phase  portraits. 


3  WHAT  IS  THE  KNOWLEDGE? 

3.1  Terminology 


The  purpose  of  this  section  is  to  introduce  some  concepts  and  definitions  from  I)y 
namical  Systems  Theory  [10].  A  dynamical  system  consists  of  two  parts:  (1)  the 


Figure  3:  Flow-chart  which  describes  the  process  of  experimenting  with  a  dynamical 
system. 

system  state,  and  (2)  the  evolution  law.  The  system  state  at  any  time  to  is  a  mini¬ 
mum  set  of  values  of  variables  {aq, . . . ,  xnj  which,  along  with  the  input  to  the  system 
for  t  >  t0,  is  sufficient  to  determine  the  behavior  of  the  system  for  all  time  t  >  t0.  The 
variables  which  define  the  system  state  are  called  state  variables.  The  conceptual 
n-dimensional  space  with  the  n  state  variables  as  basis  vectors  is  called  the  phase 
space.  A  state  vector  is  a  set  of  state  variables  considered  as  a  vector  in  the  phase 
space.  As  the  system  evolves  with  time,  the  state  vector  traces  out  a  path  in  the 
phase  space;  the  path  is  called  an  orbit  or  a  trajectory.  Finally,  a  phase  portrait 
is  a  partition  of  the  phase  space  into  orbits. 

The  evolution  law  determines  how  the  state  vector  evolves  with  time.  In  a  finite 
dimensional  discrete  time  system,  the  evolution  law  is  given  by  difference  equations. 
The  difference  equation  is  specified  by  a  function  /  :  X  — ►  X  where  X  is  the  phase 
space  of  the  discrete  system.  The  function  f  which  defines  a  discrete  dynamical  system 
is  called  a  mapping,  or  a  map,  for  short.  The  multipliers  of  the  map  f  are  the 
eigenvalues  of  the  Jacobian  of  f.  An  area-preserving  map  is  a  map  whose  Jacobian 
has  a  unit  determinant. 


The  set  of  iterates  of  f  -  x,  /(x),  /2(x),  /3(x), /n(x)  -  as  n  becomes  large  is 
called  the  orbit  of  x  relative  to  f;  it  captures  the  history  of  x  as  f  is  iterated. 

Two  types  of  point  have  the  simplest  histories  -  fixed  point,  and  periodic  point. 
The  point  x  is  a  fixed  point  of  f  if  f(x)  =  x.  A  fixed  point  x  is  called  stable,  or 
elliptic,  if  all  the  multipliers  of  f  at  x  lie  on  the  unit  circle;  it  is  called  unstable,  or 
hyperbolic,  otherwise.  The  point  x  is  a  periodic  point  of  period  n  if  fn(x)  =  x. 
The  least  positive  n  for  which  fn(x)  =  x  is  called  the  period  of  x.  The  set  of  all 
iterates  of  a  periodic  point  forms  a  periodic  orbit. 


3.2  Analytical  Techniques 

The  first  step  in  analyzing  any  dynamical  system  is  to  find  its  equilibria  and  their 
stability.  For  a  second  order  discrete  system,  this  amounts  to  solving  a  pair  of  non¬ 
linear  algebraic  equations.  To  take  an  example,  the  fixed  points  of  quadratic  map  are 
the  solutions  to  the  following  pair  of  algebraic  equations: 

x  —  icoso  —  (y  —  x2)  sin o: 

y  =  xsina  -f  (y  —  x2)cosa  (2) 

There  are  two  solutions:  x  =  y  =  0  and  x  =  2tan(y),y  =  2tan2(~-). 

Once  the  fixed  points  are  found,  the  trace  of  the  Jacobian  of  the  map  at  each  fixed 
point  is  computed  to  determine  its  linear  stability.  The  stability  condition  is  that  the 
absolute  value  of  the  trace  of  the  Jacobian  is  less  than  or  equal  to  two  [7].  For  the 
quadratic  map,  the  Jacobian  is  given  by: 

(cos a  +  2xsina  —sin a  \ 
sin  a  —  2x  cos  a  cos  a  / 

At  (0,0),  the  trace  of  the  Jacobian  is  2  cos  a.  Since  |  2  cos  a  |  <  2,  (0,0)  is  always 
stable.  Similar  stability  analysis  shows  the  other  fixed  point  is  always  unstable. 


3.3  Qualitative  Dynamics  and  their  Geometry 


3.3.1  Types  of  Orbit 


In  this  section,  I  give  a  brief  outline  of  the  most  important  types  of  orbits  in  area- 
preserving  maps.  One  important  physical  situation  that  caves  rise  to  area-preserving 
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map  is  conservative  systems  of  two  degrees  of  freedom.  Restricted  to  constant  en¬ 
ergy  surface,  the  4-dimensional  phase  space  becomes  3-dimensional.  By  the  Liouville 
Theorem  on  integrable  systems  (1,  Chapter  10],  this  3-dimensional  subspace  is  dif- 
feomorphic  to  a  family  of  concentric  2-dimensional  tori  lying  inside  one  another  (see 
figure  4).  A  torus  is  called  invariant  if  an  orbit  starting  at  a  point  on  the  surface  of 
the  torus  stays  on  it  forever. 


Figure  4:  3-dimensional  phase  space  filled  with  2-dimensional  concentric  tori.  A 
cross-section  E  is  drawn. 

Qualitatively,  three  special  types  of  motion  are  important:  (1)  periodic  motion 
lying  on  an  invariant  torus,  (2)  quasiperiodic  motion  that  densely  covers  an  in¬ 
variant  torus,  and  (3)  chaotic  motion  that  wanders  in  a  3-dimensional  volume  of 
th?  phase  space. 

We  construct  an  area-preserving  map  by  taking  a  2-dimensional  cross-section  £  of 
the  family  of  tori.  An  orbit  beginning  in  E  returns  to  it  after  making  a  circuit  around 
the  torus.  The  successive  intersections  of  the  orbit  with  E  defines  ;i  map  from  E  onto 
itself.  The  map  is  area- preserving  because  the  3-dimensional  flow  is  conservative  [2, 
Appendix  31]. 

Now,  how  do  the  periodic  and  quasiperiodic  orbits  appear  on  E?  When  the  mot  ion 
is  periodic,  the  orbit  closes  after  a  finite  number  of  intersections  with  E.  resulting  in 
a  finite  number  of  points  on  E.  If  the  motion  is  quasiperiodic,  it  will  intersect  E 
an  infinite  number  of  times,  covering  a  closed  curve  densely.  Finally,  if  the  motion 
is  chaotic,  the  orbit  will  intersect  E  erratically  at  many  different  places,  eventually 
covering  a  finite  area. 

Including  the  case  of  unbounded  motion  for  the  sake  of  completeness,  we  have  the 
following  correspondence  between  special  types  of  orbit  in  the  3-dimensional  phase 
space  and  their  geometric  manifestations  on  the  corresponding  2-dimensional  phase 
space  of  the  map: 


periodic  motion  <==>  periodic  points 

quasiperiodic  motion  *$==>  invariant  curve 
chaotic  motion  <=>  chaotic  region 

unbounded  motion  escape  points 

3.3.2  Bifurcations:  Qualitative  changes  in  the  phase  portrait 

Two  phase  portraits  are  qualitatively  equivalent  if  there  exist'  a  homeomorphism 
between  them  which  preserves  fixed  points,  periodic  points,  invariant  curves,  and 
their  stability.  Bifurcation  is  said  to  occur  when  the  dynamical  system  goe3  through 
a  qualitative  change  in  its  phase  portrait  as  the  control  parameter  is  varied.  I  will 
focus  on  one  important  type  of  bifurcation:  appearance  and  disappearance  of  periodic 
orbits. 

I  begin  with  [12]  which  gives  a  complete  classification  of  the  generic  bifurcations  of 
periodic  points  for  one-parameter  area-preserving  maps.  Only  bifurcations  which  are 
generic  are  classified  because  the  pathological  cases  arise  rarely,  and  are  much  harder 
to  deal  with.  Let  A  and  j  be  the  multipliers  of  a  periodic  point  of  an  area-preserving 
map.  Meyer  has  shown  that  generic  bifurcations  occur  when  A  =  nth  root  of  unity 
where  n  =  1,  2,  3,  4,  and  >  5.  Let  us  consider  each  of  these  cases  (see  figure  5): 

Case  1.  A  =  1  Extremal  Bifurcation 

The  (discrete)  flow  line  forms  a  cusp  at  the  bifurcation  point  and  as  the  multiplier 
passes  through  +1,  there  is  a  sudden  birth  of  a  pair  of  elliptic  and  hyperbolic  closed 
orbits.  Thus,  a  pair  of  fixed  points  are  created  for  the  map.  This  is  also  known  as 
creation,  saddle-node  bifurcation,  or  tangent  bifurcation. 

Case  2.  A  —  —  1  Transitional  Bifurcation 

As  the  multiplier  passes  through  -1,  the  stable  elliptic  periodic  point  loses  its 
stability,  i.e.,  turns  hyperbolic,  and  a  stable  cycle  of  twice  the  period  appears.  Other 
names  for  this  phenomenon  are:  period  doubling,  flip  bifurcation,  subtle  division. 

Case  3.  A  =  e~  Phantom  3-kiss 


The  region  of  stability  of  the  elliptic  fixed  point  shrinks  to  zero  as  the  hyperbolic 
points  of  an  unstable  period-3  cycle  “kiss”  at  the  origin.  After  the  “kiss”,  the  fixed 
point  turns  elliptic  again,  and  a  new  unstable  period-3  cycle  is  emitted.  Note  the 
change  in  orientation  of  the  triangular  region  around  the  elliptic  point.  The  phantom 
3-kiss  is  often  preceded  by  extremal  bifurcations  in  a  region  a  bit  further  away  from 
the  original  elliptic  fixed  point,  resulting  in  the  formation  of  a  pair  of  elliptic  and 
hyperbolic  period-3  points. 

Case  4.  A  =  e2?1  Phantom  4-kiss 


When  the  multiplier  A  is  a  4-th  root  of  unity,  it  can  undergo  two  types  of  bifurca¬ 
tion.  The  first  case,  phantom  4-kiss,  is  similar  to  that  of  case  3  except  that  we  now 
have  hyperbolic  points  of  a  period-4  orbit  “kissing”.  Tangent  bifurcation  also  occurs 
at  a  distance  from  the  center  elliptic  fixed  point,  resulting  in  a  an  unstable  period-4 
cycle.  The  second  case  belongs  to  the  fifth  type  to  be  described  next. 

2frt 

Case  5.  A  ==  e~  where  p  >  5  Emission 

As  the  multiplier  passes  through  a  p-th  root  of  unity  where  p  >  5,  there  is  an 
emission  of  elliptic  and  hyperbolic  sub-harmonics  of  k  0f  the  orbital  frequency  of  the 
original  elliptic  orbit  (shown  in  the  figure  for  the  case  p  =  5).  The  elliptic  points  of 
the  stable  period-p  cycle  alternate  with  the  hyperbolic  points  of  the  unstable  one  near 
the  vertices  of  a  regular  p-gon  with  center  at  origin.  Thus,  emission  is  responsible  for 
the  creation  of  “island  chains”  in  the  phase  portraits  of  Figure  1.  Meyer  called  this 
type  of  bifurcation  the  generic  p-bifurcation. 


4  APPROACH  TO  THE  CONTROL  PROBLEM 

4.1  How  to  start  the  numerical  experiment? 

Elliptic  fixed  points  are  good  places  to  start.  We  expect  that  the  orbits  near  an 
elliptic  fixed  point,  where  the  linear  terms  of  the  map  dominate,  will  be  mostly 
invariant  curves.  We  then  search  radially  outward  until  we  encounter  island  chains, 
and  eventually  chaotic  regions. 


AMV.V.V 


Figure  5:  Five  Generic  types  of  Bifurcation  Geometry.  A  periodic  point  bifurcates 
whenever  its  multipliers  pass  through  an  n-th  root  of  unity.  The  bifurcation  diagram 
shows  the  loci  of  fixed  points  as  the  parameter  t  is  varied. 

4.2  How  to  decide  what  experiment  to  try  next? 

Knowing  the  generic  bifurcation  patterns  is  valuable  for  controlling  numerical  exper¬ 
iments.  To  begin  with,  it  is  difficult  to  locate  the  value  of  the  control  parameter 
at  which  bifurcation  occurs:  it  is  of  probability  almost  zero  that  a  randomly  chosen 
point  in  the  (x,  y ,  a)  space  will  be  the  bifurcation  point.  But  the  pattern  of  flow  near 
a  periodic  point  just  before  and  after  the  bifurcation  occurs  in  a  finite  range  of  the 
control  parameter;  hence  the  pattern  is  easier  to  detect  during  the  experiments. 

Once  a  given  flow  pattern  is  found  to  match  some  parts  in  our  library  of  bifurcation 
geometries,  it  will  give  us  strong  evidence  that  the  corresponding  bifurcation  exists, 
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and  we  should  be  able  to  locate  the  rest  of  the  flow  patterns  as  given  by  the  generic 
bifurcation.  The  pre-stored  knowledge  about  these  bifurcations  gives  us  the  complete 
information  about  what  geometric  objects,  and  approximately  where  in  the  control 
parameter  space  to  look  for. 

To  take  an  example,  consider  the  phantom  3-kiss  seen  in  figure  Id.  The  local  flow 
pattern  around  the  fixed  point  matches  that  in  figure  5.  According  to  the  bifurcation 
pattern,  the  regular  region  around  the  stable  fixed  point  will  shrink  in  size,  becoming 
an  unstable  fixed  point;  eventually,  a  new  stable  fixed  point  is  born.  So,  we  should 
expect  to  see  figure  If  at  some  a  slightly  greater  than  two. 

4.3  How  to  decide  when  to  terminate  the  experiments? 

Besides  imposing  a  strong  constraint  on  what  can  be  expected  to  happen  in  the  phase 
portrait,  the  generic  bifurcations  also  provide  an  answer  to  the  problem  of  termina¬ 
tion:  a  simulation  experiment  is  incomplete  unless  all  the  major  qualitative  features 
in  the  phase  portrait  can  be  explained  by  this  finite  list  of  local  generic  bifurcations. 
An  example  is  the  change  of  stability  of  a  fixed  point.  Suppose  we  have  numerically 
located  the  3-island  chain  and  the  center  elliptic  point  at  some  a  =  a0  as  in  figure 
Id.  We  know  that  the  family  of  phase  portraits  is  yet  incomplete  because  we  expect 
a  phantom  3-kiss  bifurcation  to  occur.  In  particular,  we  need  to  try  at  least  two  more 
experiments  to  obtain  two  phase  portraits:  first,  at  a  =  aj  when  the  triangular  region 
is  flipped,  and  second,  at  a2  6  (a0,  Qi)  when  the  region  becomes  vanishingly  small, 
indicating  instability  of  the  fixed  point. 

5  SOLVING  THE  INTERPRETATION  PROBLEM 

The  Interpretation  Problem  consists  of  the  following  sub-problems: 

1.  Orbit  Type.  How  can  one  recognize  the  orbit  type  -  a  O-dimensional  finite  point 
set  whose  elements  are  encountered  repeatedly,  a  1 -dimensional  smooth  curve, 
or  a  2-dimensional  region  -  of  a  set  of  iterates? 

2.  Clustering.  How  can  one  determine  the  number  of  islands  in  an  island  chain? 
This  number  gives  the  period  of  the  enclosed  periodic  point. 


3.  Area  and  Centroid.  How  can  one  estimate  the  centroid  and  area  enclosed  by 
the  curve?  The  centroid  is  a  good  approximation  of  the  location  of  the  enclosed 
periodic  point.  The  area  gives  a  measure  of  saliency  of  the  island  chain. 

4.  Shape.  How  can  one  recognize  the  shape  of  a  curve?  For  example,  is  it  a  3-sided 
figure  resembling  a  triangle? 

In  the  following,  I  will  show  how  these  four  problems  can  be  solved  by  applying 
techniques  from  computational  geometry  and  computer  vision. 

5.1  Euclidean  Minimal  Spanning  Tree 

Given  n  points  in  the  Euclidean  plane,  a  Euclidean  minimal  spanning  tree  (EM ST)  is 
a  tree  that  interconnects  all  the  points  with  minimal  total  edge  length.  It  turns  out 
that  EMST  provides  us  with  the  fundamental  data  structure  to  solve  the  first  three 
sub-problems. 

Figure  6  shows  the  steps  of  processing.  First,  the  program  computes  a  EMST  from 
the  input  point  set.  The  current  implementation  uses  the  Prim-Dijkstra  algorithm 
which  has  a  run  time  complexity  of  0(n2)  where  n  is  the  number  of  points  [3].  It  takes 
about  100  seconds  on  the  Symbolics  to  compute  a  EMST  for  a  set  of  512  points.  A 
faster  algorithm  due  to  Shamos  [13]  which  makes  use  of  the  Voronoi  diagram  is  being 
implemented.  The  new  algorithm  can  achieve  the  optimal  lower  bound  of  O (n  logn) 
time.  Figure  7  shows  the  EMSTs  constructed  for  three  sample  point  sets. 

Second,  the  program  detects  clusters  in  the  EMST  by  looking  for  edges  in  the  tree 
that  are  significantly  longer  than  nearby  edges.  Such  edges  are  called  inconsistent 
[17].  The  following  criterion  suggested  by  Zahn  is  used  to  determine  if  an  edge  is 
inconsistent.  An  edge  is  inconsistent  if  (1)  its  length  is  more  than  two  times  the 
average  of  nearby  edges,  and  (2)  its  length  is  more  than  two  standard  deviations 
larger  than  the  average  of  the  lengths  of  nearby  edges  (where  the  statistics  are  taken 
from  the  set  of  all  nearby  edges).  A  point  P  is  nearby  point  Q  if  P  is  connected  to 
Q  in  the  EMST  containing  two  or  fewer  edges.  Inconsistent  edges  are  then  deleted. 
This  breaks  up  the  EMST  into  connected  sub-components  which  are  collected  by  a 
depth-first  tree  walk.  Figure  8  shows  the  result  of  the  clustering  procedure  on  the 
input  point  set  of  figure  7b.  The  procedure  finds  three  inconsistent  edges.  For  the 
other  two  point  sets,  no  inconsistent  edge  is  detected.  The  output  of  the  clustering 


INPUT:  A  POINT  SET 


Figure  6:  Main  processing  steps  for  solving  the  interpretation  problem 


procedure  is  a  collection  of  connected  sub-components  each  of  which  represents  a 
sub- tree  of  the  original  EMST. 


Third,  for  each  sub-tree  of  the  EMST,  the  program  examines  the  degree  of  each 
of  its  nodes,  where  the  degree  of  a  node  is  the  number  of  nodes  connected  to  it  in 
the  sub-tree.  For  a  smooth  curve,  the  EMST  consists  of  two  terminal  nodes  of  degree 
one;  the  rest,  degree  two.  For  a  point  set  that  fills  an  area,  its  corresponding  EMST 
consists  of  many  nodes  having  degree  three  or  higher. 


Finally,  to  compute  the  area  and  centroid  of  the  region  bounded  by  a  curve, 
the  program  generates  an  ordered  sequence  of  points  from  the  EMST,  and  spline- 
interpolates  the  sequence  to  obtain  a  smooth  curve.  The  smooth  curve  is  encoded 
using  chain  coding  [6].  Straightforward  algorithms  are  then  applied  to  compute  the 
area  and  centroid. 
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Figure  8:  Result  of  clustering.  Three  inconsistent  edges  are  detected  in  the  point  set 
of  figure  7(b). 


5.2  Scale  Space  Image 


The  shape  of  a  curve,  as  we  have  seen  in  section  4,  is  an  important  clue  to  determine 
what  type  of  bifurcation  is  occurring.  The  number  of  sides,  the  locations  and  types 
of  extremal  curvature  -  these  are  two  of  the  most  important  geometric  properties  to 
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recognize. 


The  recognition  algorithm  is  based  on  the  method  of  scale  space  image  introduced 
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by  Witkin  [16].  First,  a  curve  is  parameterized  by  C(s)  =  (x(s),y(s))  where  s  is 
the  arc  length  along  the  curve.  The  two  functions  x(s)  and  y(s)  are  periodic  for  a 
closed  curve,  and  they  can  be  computed  from  the  chain  code  representation.  Second, 
x(s)  and  y(s)  are  smoothed  by  the  Gaussian  and  its  first  two  derivatives  at  multiple 
spatial  scales.  The  scales  are  separated  by  one  octave:  4,  5,  7,  11,  15,  and  22,  which 
are  results  of  successively  multiplying  4  by  \/2  and  rounding.  Third,  the  curvature 
function  <c(s)  is  computed  by  the  formula: 

_  xy  -  yx 

K  ~  (i2  +  y 2)i 

Finally,  the  zero-crossings  of  #c(s),  and  the  signs  of  «(s)  are  computed  to  determine 
the  locations  and  type  of  the  extrema. 
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Figure  9:  Type  a  nd  locations  of  extrema  of  curvature.  Spatial  scales:  4,  5,  7,  11,  15, 
and  22.  •  =  maximum  curvature;  x  =  minimum  curvature. 

An  example  of  the  operation  of  the  recognition  algorithm  is  shown  in  Figure  9. 
The  input  is  an  orbit  segment  of  256  points  which  seems  to  fill  a  3-sided  curve.  The 
positions  of  the  extrema  of  curvature  found  at  various  smoothing  scales  are  marked: 
•  for  maximum  curvature;  x  for  minimum  curvature. 


6  Summary 


In  this  paper,  I  have  studied  the  task  of  qualitative  analysis  of  nonlinear  area¬ 
preserving  map  by  numerical  experiments.  I  have  also  described  how  to  approach 
the  two  major  problems  in  automating  the  experimenting  process:  (1)  experiment 
control,  and  (2)  result  interpretation.  The  basic  idea  is  that  knowledge  of  qualitative 
dynamics  and  bifurcations  provides  a  strong  constraint  on  the  type  of  behavior  pos- 
.  sible.  Finally,  I  have  described  a  program  which  solves  the  interpretation  problem  by 
using  techniques  from  computational  geometry  and  computer  vision. 
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